rm(list = ls())
library(foreign)
library(sandwich)
library(CBPS)
library(MASS)
library(rms)
library(xtable)
library(plyr)
library(fmsb)
library(stargazer)
library(stats)
library(mice)
library(Amelia)
set.seed(312349)

#Set working directory and load data
load('data.rData')

#Useful functions
checkBinaryTrait <- function(v, naVal="NA") {
	if (!is.numeric(v)) stop("Only numeric vectors are accepted.")
	vSet = unique(v)
	if (!missing(naVal)) vSet[vSet == naVal] = NA
	vSet = vSet[!is.na(vSet)]

	if (any(as.integer(vSet) != vSet)) FALSE
	else if (length(vSet) > 2) FALSE
	else TRUE
}
deci <- function(x, k, center = TRUE) {
	a <- format(round(x, k), nsmall=k)
	return(a)
}
decim <- function(x, k, center = TRUE) {
	a <- format(round(x, k), nsmall=k)
	if (center == TRUE & x >= 0.0){
		a <- paste('**', a, sep = '')
	}
	return(a)
}
testStatVector <- function(dat, balanceVars){
	returnVector <- NULL
	for (var in balanceVars){
		test.stat <- mean(dat[dat$cp == 1, var], na.rm = T) - mean(dat[dat$cp == 0, var], na.rm = T)
		returnVector <- c(returnVector, test.stat)
	}
	names(returnVector) <- balanceVars
	return(returnVector)
}

boot_lm <- function(dat, formul, cluster_var = 'ID', strata = NULL, R = 500){
	sample_coefs = NULL
	the_clusters <- as.character(unique(dat[[cluster_var]]))
	ref.model <- lm(formul, dat)
	coef_labels <- names(ref.model$coefficients)
	for (i in 1:R){
		boot_data <- data.frame()
		states.selected <- sample(x = the_clusters, size = length(the_clusters), replace = TRUE)
		bb <- table(states.selected)
		for(j in 1:max(bb)){
			selected.rows <- dat[[cluster_var]] %in% names(bb[bb %in% j])
			add_row <- dat[selected.rows, ]
			for(k in 1:j){
				boot_data <- rbind(boot_data, add_row)
			}
		}
		newrow <- rep(NA, length(coef_labels))
		names(newrow) <- coef_labels

		if(!is.null(strata)){
			boot_data$lm_weights <- NULL
			strat_levels <- levels(as.factor(boot_data[[strata]]))
			for (strat_level in strat_levels){
				this_strat <- boot_data[boot_data[[strata]] == strat_level,]
				pscore <- mean(this_strat$cp, na.rm = T)
				boot_data$lm_weights[boot_data[[strata]] == strat_level & boot_data$cp == 1] <- 1/pscore
				boot_data$lm_weights[boot_data[[strata]] == strat_level & boot_data$cp == 0] <- 1/(1-pscore)
			}
			#lm.out <- lm(formula = update(formul, . ~ . - factor(nAttempts)), data = boot_data, weights = lm_weights)
			lm.out <- lm(formula = formul, data = boot_data, weights = lm_weights)
			}else{lm.out <- lm(formula = formul, data = boot_data)}
		the.coefs <- summary(lm.out)$coefficients[, 1]
		for (thevar in names(the.coefs)){
			newrow[names(newrow) == thevar] <- the.coefs[names(the.coefs) == thevar]
		}
		sample_coefs <- rbind(sample_coefs, newrow)
	}
	rownames(sample_coefs) <- 1:R
	return(sample_coefs)
}
lm_mi <- function(dat_t, controls = NULL, reweigh = NULL, nAttempts_cntrl = TRUE, ctrend = F){
	coefs <- NULL
	ses <- NULL
	if (!is.null(controls)){
		reg.formul <- formula(paste('DV_diff ~ cp + factor(nAttempts)', ' + ', paste(controls, '.pt1', sep = '', collapse = ' + '), sep = ''))
	}else{
		reg.formul <- DV_diff ~ cp + factor(nAttempts)
	}
	if (nAttempts_cntrl == FALSE){
		reg.formul <-  update(reg.formul, . ~. -factor(nAttempts))
	}
	if (is.null(reweigh)){
		for (i in levels(dat_t$.imp)) {
		    in.imp <- dat_t$.imp == i
		    dt <- dat_t[in.imp, ]
		    did.out <- lm(reg.formul, data = dt, weights = NULL)
			new_coef <- summary(did.out)$coef[,1]
			new_se <- summary(did.out)$coef[,2]
			coefs <- rbind(coefs, new_coef)
			ses <- rbind(ses, new_se)
		}
	}else{
		for (i in levels(dat_t$.imp)) {
		    in.imp <- dat_t$.imp == i
		    dt <- dat_t[in.imp, ]
		    did.out <- lm(reg.formul, data = dt, weights = dt[[reweigh]])
			new_coef <- summary(did.out)$coef[,1]
			new_se <- summary(did.out)$coef[,2]
			coefs <- rbind(coefs, new_coef)
			ses <- rbind(ses, new_se)
		}
	}
	coefs_combined <- as.vector(mi.meld(coefs, ses)$q.mi)
	ses_combined <- as.vector(mi.meld(coefs, ses)$se.mi)
	zvals_t <- -as.vector(abs(coefs_combined/ses_combined))
	zTests <- 2*pnorm(zvals_t)

	a <- cbind(coefs_combined, ses_combined, 2*pnorm(zvals_t))
	rownames(a) <- colnames(mi.meld(coefs, ses)$q.mi)
	colnames(a) <- c('Coefs', 'Ses', 'PVal')
	return(list('coefs' = a, 'N' = nrow(dt), 'Nsucc' = sum(dt$cp), model = did.out))
}
vcovCluster <- function(model, cluster){
	if(nrow(model.matrix(model))!=length(cluster)){
	  stop("Cluster variable and model have different N")
	}
	M <- length(unique(cluster))
	N <- length(cluster)
	K <- model$rank
	dfc <- (M/(M - 1)) * ((N - 1)/(N - K ))
	uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
	rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
	return(rcse.cov)
}
panel_reg <- function(data.in, formul, ...){
	arguments <- list(...)
	if ('weights' %in% names(arguments)){
		if (!is.null(arguments[['weights']])){
			weights.in <- names(arguments)
			weights.lab <- unlist(arguments)
			data.in$weights <- data.in[[weights.lab]]
			did.out <- lm(formula = as.formula(formul), data = data.in, y = TRUE, weights = weights)
		}else{
			did.out <- lm(formula = as.formula(formul), data = data.in, y = TRUE, weights = NULL)
		}
	}else{
		did.out <- lm(formula = as.formula(formul), data = data.in, y = TRUE, weights = NULL)
	}
	#lm(as.formula(formul), data = data.in, y = TRUE, weights = data.in[[cus_weights]])
	robust_vcov <- vcovCluster(model = did.out, cluster = data.in$ccode)
	a <- coeftest(did.out, robust_vcov, df.residual(did.out))
	return(list(model = did.out, coefs = a))
}

#-------------------------------
#Section A2
#-------------------------------
#Non-Parametric combination test
balanceVars <- c(
	'e_migdppclnL1', #GDP per capita
	'e_migdpgroL1', #GDP growth
	'ln.rexportspgdpL1', #Exports as share of GDP (log)
	'oil_binaryL1', #Oil and gas production
	'ethfracL1', #Ethnic fractionalization
	'ln.tpopL1', #Total population (logged)
	'tpopgroL1', #Population growth
	'cheibubDemoL1', #Democracy (CGV)
	'cheibubPresidL1', #Presidentialism (CGV)
	'cheibubMilL1',  #Military regime (CGV)
	'gwf_partyL1', #Party regime (GWF)
	'gwf_monarchL1', #Monarchy (GWF)
	'gwf_personalL1', #Personalist regime (GWF)
	'cheibubLegisL1', #Legislature (CGV)
	'ln.milperpcL1', #Per capita military personnel (log)
	'ln.milperL1', #Total military personnel (log)
	'ln.rCOWmilexL1', #Total military expenditures (log)
	'ln.CowmilexpgdpL1', #Military expenditures, share of GDP (log)
	'ln.rCOWmilexpcL1', #Per capita military expenditures (log)
	'ln.rCOWmilexpsoldL1', #Per soldier military expenditures (log)
	'cbd.changeCowMilexL1', #Change in military expenditures (cubic root)
	'effectivenumberL1', #Counterbalancing,  Pilster and Bohmelt
	'paramilitary.bscL1', #Counterbalancing, Belkin and Schofer
	'counterbal.pthL1', #Counterbalancing, Powell and Thyne
	'counterbalancingL1', #Counterbalancing, De Bruin
	'warL1', #Interstate war
	'CW', #Cold War
	'intra_warL1', #Intra-state war
	'prioL1', #UCDP Prio
	'OppNonViolBanksL1', #Unrest, non violent (Banks)
	'ncp.bsc', #Number of past successful coups
	'ncpFl.bsc', #Number of past failed coups
	'purgesL1', #Purges
	'latentmeanL1'# State repression
	 )

monte_carlo_dat <- dat[, c(balanceVars)]
prob_treatment <- sum(dat$cp)/nrow(dat)
simulations <- 1:1000
testStatSims <- NULL
for (simulation in simulations){
	monte_carlo_dat <- dat
	treated <- sample(1:nrow(dat), round(nrow(dat)*prob_treatment))
	monte_carlo_dat$cp <- 0
	monte_carlo_dat$cp[treated] <- 1
	testStatSims <- rbind(testStatSims, testStatVector(monte_carlo_dat, balanceVars))
}

Tobs <- testStatVector(dat, balanceVars)
Lobs <- NULL
Lobs_star <- NULL
for (var in balanceVars){
	Lobs_t <- 1-sum(abs(testStatSims[, var]) < abs(Tobs[var]))/nrow(testStatSims)
	Lobs <- c(Lobs, Lobs_t)

	Lobs_tar_t <- apply(as.matrix(testStatSims[, var]), 1, function(x) sum(x < testStatSims[, var])/nrow(testStatSims))
	Lobs_star <- rbind(Lobs_star, Lobs_tar_t)
}
names(Lobs) <- balanceVars
rownames(Lobs_star) <- balanceVars
TobsGlobal <- -2*sum(log(Lobs))
TobsGlobal_star = apply(Lobs_star, 2, function(x) -2*sum(log(x)))
NPC_test <- sum(TobsGlobal_star>TobsGlobal)/length(TobsGlobal_star)

#-------------------------------
#TABLE A1
#-------------------------------
library(lmtest)

controls <- c(
			'cheibubMil',
			'ythblgap',
			'cheibubDemo',
	 		'ln.rCOWmilexpc',
			'CW',
			'counterbalancing',
			'ncp.bsc',
			'ncpFl.bsc'
			)
match.on <- c(
	 'e_migdppcln',
	 'cheibubLegis',
	 'ln.milperpc',
	 'counterbalancing',
	 'oil_binary',
	 'purges',
	 'ncp.bsc',
	 'ncpFl.bsc',
	 'OppNonViolBanks',
	 'ethfrac'
	 )

diff_dat$nAttempts[diff_dat$nAttempts >= 3] <- 3
diff_dat$ID <- 1:nrow(diff_dat)
diff_temp1 <- subset(diff_dat, nAttempts ==1)
diff_temp2 <- subset(diff_dat, nAttempts ==2)
diff_temp3 <- subset(diff_dat, nAttempts ==3)
pscore1 <- mean(diff_temp1$cp)
pscore2 <- mean(diff_temp2$cp)
pscore3 <- mean(diff_temp3$cp)

diff_dat$lm_weights[diff_dat$nAttempts == 1 & diff_dat$cp == 1] <- 1/pscore1
diff_dat$lm_weights[diff_dat$nAttempts == 1 & diff_dat$cp == 0] <- 1/(1-pscore1)
diff_dat$lm_weights[diff_dat$nAttempts == 2 & diff_dat$cp == 1] <- 1/pscore2
diff_dat$lm_weights[diff_dat$nAttempts == 2 & diff_dat$cp == 0] <- 1/(1-pscore2)
diff_dat$lm_weights[diff_dat$nAttempts == 3 & diff_dat$cp == 1] <- 1/pscore3
diff_dat$lm_weights[diff_dat$nAttempts == 3 & diff_dat$cp == 0] <- 1/(1-pscore3)

diff_dat$ID <- as.character(1:nrow(diff_dat))

diff_dat_mi$ID <- 1:nrow(diff_dat_mi)
diff_dat_mi$ccode <- factor(as.character(diff_dat_mi$ccode))
diff_dat_mi$w <- NA
match.formul <- formula(paste('cp ~ ', paste(paste(match.on, '.pt1', sep = '', collapse = ' + '), sep = ''), sep = ''))
for (i in levels(diff_dat_mi$.imp)) {
	in.imp <- diff_dat_mi$.imp == i
	m.out <- CBPS(match.formul, standardize = F, ATT = 0, data = diff_dat_mi[in.imp,])
	diff_dat_mi$w[in.imp] <- m.out$weights
	diff_dat_mi$ID[in.imp] <- diff_dat$ID
}

#------------------------------------------
#Baseline
formul <- DV_diff ~ cp + factor(nAttempts)
did.out_boot <- boot_lm(dat = diff_dat, formul = formul, cluster_var = 'ID', strata = NULL)
did.out.nostrat <- lm(formul, data = diff_dat, weights = NULL)
coefs.out.nostrat <- coef(did.out.nostrat)
SEs.out.nostrat <- apply(did.out_boot, 2, sd)
pvalues.out.nostrat <- 2*pnorm(-abs(coefs.out.nostrat/SEs.out.nostrat))

estims.cov <- lm_mi(dat_t = diff_dat_mi, controls = controls, reweigh = 'w', nAttempts_cntrl = TRUE)
coefs.cov.out.nostrat <- estims.cov$coefs[,1]
SEs.cov.out.nostrat <- estims.cov$coefs[,2]
pvalues.cov.out.nostrat <- estims.cov$coefs[,3]
N.cov.out.nostrat <- estims.cov$N
Nsucc.cov.out.nostrat <- estims.cov$Nsucc
did.out.cov.nostrat <- estims.cov$model

#Matching
#------------------------------------------
estims.match <- lm_mi(dat_t = diff_dat_mi, controls = NULL, reweigh = 'w')
coefs.match.out.nostrat <- estims.match$coefs[,1]
SEs.match.out.nostrat <- estims.match$coefs[,2]
pvalues.match.out.nostrat <- estims.match$coefs[,3]
N.match.out.nostrat <- estims.match$N
Nsucc.match.out.nostrat <- estims.match$Nsucc
did.out.match.nostrat <- estims.match$model

estims.match.cov <- lm_mi(dat_t = diff_dat_mi, controls = controls, reweigh = 'w')
coefs.match.cov.out.nostrat <- estims.match.cov$coefs[,1]
SEs.match.cov.out.nostrat <- estims.match.cov$coefs[,2]
pvalues.match.cov.out.nostrat <- estims.match.cov$coefs[,3]
N.match.cov.out.nostrat <- estims.match.cov$N
Nsucc.match.cov.out.nostrat <- estims.match.cov$Nsucc
did.out.match.cov.nostrat <- estims.match.cov$model

#Trends
#------------------------------------------
formul_trend <- DV_diff ~ cp + factor(nAttempts) + trend.pt1
did.out_boot_trend <- boot_lm(dat = diff_dat, formul = formul_trend, cluster_var = 'ID', strata = NULL, R = 500)
did.out_trend <- lm(formul_trend, data = diff_dat, weights = NULL)
coefs.trends.out.nostrat <- coef(did.out_trend)
SEs.trends.out.nostrat <- apply(did.out_boot_trend, 2, sd)
pvalues.trends.out.nostrat <- 2*pnorm(-abs(coefs.trends.out.nostrat/SEs.trends.out.nostrat))
N.trend.out.nostrat <- nrow(diff_dat)
Nsucc.cov.trend.out.nostrat <- sum(diff_dat$cp)
did.trend.out.nostrat <- did.out_trend

estims.trendsqrd <- lm_mi(dat_t = diff_dat_mi, controls = c(controls, 'trend', 'trendsqrd'), reweigh = 'w')
coefs.trendsSqrd.out.nostrat <- estims.trendsqrd$coefs[,1]
SEs.trendsSqrd.out.nostrat <- estims.trendsqrd$coefs[,2]
pvalues.trendsSqrd.out.nostrat <- estims.trendsqrd$coefs[,3]
N.trendssqr.out.nostrat <- estims.trendsqrd$N
Nsucc.trendssqr.out.nostrat <- estims.trendsqrd$Nsucc
did.out.trendssqr.nostrat <- estims.trendsqrd$model

coefsList <- list(
	coefs.out.nostrat,
	coefs.cov.out.nostrat,
	coefs.match.out.nostrat,
	coefs.match.cov.out.nostrat,
	coefs.trends.out.nostrat,
	coefs.trendsSqrd.out.nostrat
	)

seList <- list(
	SEs.out.nostrat,
	SEs.cov.out.nostrat,
	SEs.match.out.nostrat,
	SEs.match.cov.out.nostrat,
	SEs.trends.out.nostrat,
	SEs.trendsSqrd.out.nostrat
	)

pList <- list(
	pvalues.out.nostrat,
	pvalues.cov.out.nostrat,
	pvalues.match.out.nostrat,
	pvalues.match.cov.out.nostrat,
	pvalues.trends.out.nostrat,
	pvalues.trendsSqrd.out.nostrat
	)

successes <- c(
	sum(diff_dat$cp),
	Nsucc.cov.out.nostrat,
	Nsucc.match.out.nostrat,
	Nsucc.match.cov.out.nostrat,
	Nsucc.cov.trend.out.nostrat,
	Nsucc.trendssqr.out.nostrat
	)

attempts <- c(
	nrow(diff_dat),
	N.cov.out.nostrat,
	N.match.out.nostrat,
	N.match.cov.out.nostrat,
	N.trend.out.nostrat,
	N.trendssqr.out.nostrat
	)

TableA1 <- stargazer(
	did.out.nostrat,
	did.out.cov.nostrat,
	did.out.match.nostrat,
	did.out.match.cov.nostrat,
	did.trend.out.nostrat,
	did.out.trendssqr.nostrat,
	coef = coefsList,
	se = seList,
	p = pList,
	p.auto = FALSE,
	omit = c('factor', 'region', 'Constant'),
	dep.var.caption  = "\\emph{Dependent variable}: First Difference in State Repression",
	dep.var.label  = c('', '', '', '', '', ''),
	column.labels   = c('Baseline Models', "With CBPS Weighing", 'With Trends'),
	column.separate = c(2, 2, 2),
	dep.var.labels.include  = FALSE,
	align = TRUE,
	omit.stat = c('rsq', 'f', 'ser', 'adj.rsq', 'n'),
	add.lines = list(
		c("CBPS Weighting", '$No$', '$No$', '$Yes$', '$Yes$', '$No$', '$No$'),#, '$No$', '$No$', '$Yes$', '$Yes$'),
		c("Observations", attempts),
		c("Successful Coups", successes)
		),
	notes.append = FALSE,
	#label = 'did.results',
	float = FALSE
	)

#-------------------------------
#TABLE A2
#-------------------------------
diff_dat_no_consec <- subset(diff_dat, nYears == 1)
diff_dat_no_consec$ID <- 1:nrow(diff_dat_no_consec)
diff_dat_no_consec$ccode <- factor(as.character(diff_dat_no_consec$ccode))

diff_dat_no_consec_mi <- subset(diff_dat_mi, nYears == 1)
diff_dat_no_consec_mi$w <- NA
match.formul <- formula(paste('cp ~ ', paste(paste(match.on, '.pt1', sep = '', collapse = ' + '), sep = ''), sep = ''))
for (i in levels(diff_dat_no_consec_mi$.imp)) {
	in.imp <- diff_dat_no_consec_mi$.imp == i
	m.out <- CBPS(match.formul, standardize = F, ATT = 0, data = diff_dat_no_consec_mi[in.imp,])
	diff_dat_no_consec_mi$w[in.imp] <- m.out$weights
	diff_dat_no_consec_mi$ID[in.imp] <- diff_dat_no_consec$ID
}

#------------------------------------------
#Baseline
formul <- DV_diff ~ cp + factor(nAttempts)
did.out_boot <- boot_lm(dat = diff_dat_no_consec, formul, cluster_var = 'ID', strata = NULL)
did.out.noconsec <- lm(formul, data = diff_dat_no_consec, weights = NULL)
coefs.out.noconsec <- coef(did.out.noconsec)
SEs.out.noconsec <- apply(did.out_boot, 2, sd)
pvalues.out.noconsec <- 2*pnorm(-abs(coefs.out.noconsec/SEs.out.noconsec))

estims.cov <- lm_mi(dat_t = diff_dat_no_consec_mi, controls = controls, reweigh = NULL)
coefs.cov.out.noconsec <- estims.cov$coefs[,1]
SEs.cov.out.noconsec <- estims.cov$coefs[,2]
pvalues.cov.out.noconsec <- estims.cov$coefs[,3]
N.cov.out.noconsec <- estims.cov$N
Nsucc.cov.out.noconsec <- estims.cov$Nsucc
estims.cov.noconsec <- estims.cov$model

#Matching
#------------------------------------------
estims.match <- lm_mi(dat_t = diff_dat_no_consec_mi, controls = NULL, reweigh = 'w')
coefs.match.out.noconsec <- estims.match$coefs[,1]
SEs.match.out.noconsec <- estims.match$coefs[,2]
pvalues.match.out.noconsec <- estims.match$coefs[,3]
N.match.out.noconsec <- estims.match$N
Nsucc.match.out.noconsec <- estims.match$Nsucc
estims.match.noconsec <- estims.match$model

estims.match.cov <- lm_mi(dat_t = diff_dat_no_consec_mi, controls = controls, reweigh = 'w')
coefs.match.cov.out.noconsec <- estims.match.cov$coefs[,1]
SEs.match.cov.out.noconsec <- estims.match.cov$coefs[,2]
pvalues.match.cov.out.noconsec <- estims.match.cov$coefs[,3]
N.match.cov.out.noconsec <- estims.match.cov$N
Nsucc.match.cov.out.noconsec <- estims.match.cov$Nsucc
estims.match.cov.noconsec <- estims.match.cov$model

#Trends
#------------------------------------------
formul_trend <- DV_diff ~ cp + trend.pt1 +  factor(nAttempts)
did.out_boot_trend <- boot_lm(dat = diff_dat_no_consec, formul_trend, cluster_var = 'ID', strata = NULL)
did.out_trend.noconsec <- lm(formul_trend, data = diff_dat_no_consec, weights = NULL)

coefs.trends.out.noconsec <- coef(did.out_trend.noconsec)
SEs.trends.out.noconsec <- apply(did.out_boot_trend, 2, sd)
pvalues.trends.out.noconsec <- 2*pnorm(-abs(coefs.trends.out.noconsec/SEs.trends.out.noconsec))
N.trend.out.noconsec <- nrow(diff_dat)
Nsucc.cov.trend.out.noconsec <- sum(diff_dat$cp)

#------------------------------------------
estims.trendsqrd <- lm_mi(dat_t = diff_dat_no_consec_mi, controls = c(controls, 'trend', 'trendsqrd'), reweigh = NULL)
coefs.trendsSqrd.out.noconsec <- estims.trendsqrd$coefs[,1]
SEs.trendsSqrd.out.noconsec <- estims.trendsqrd$coefs[,2]
pvalues.trendsSqrd.out.noconsec <- estims.trendsqrd$coefs[,3]
N.trendssqr.out.noconsec <- estims.trendsqrd$N
Nsucc.trendssqr.out.noconsec <- estims.trendsqrd$Nsucc
estims.trendsqrd.noconsec <- estims.trendsqrd$model

coefsList <- list(
	coefs.out.noconsec,
	coefs.cov.out.noconsec,
	coefs.match.out.noconsec,
	coefs.match.cov.out.noconsec,
	coefs.trends.out.noconsec,
	coefs.trendsSqrd.out.noconsec
	)

seList <- list(
	SEs.out.noconsec,
	SEs.cov.out.noconsec,
	SEs.match.out.noconsec,
	SEs.match.cov.out.noconsec,
	SEs.trends.out.noconsec,
	SEs.trendsSqrd.out.noconsec
	)

pList <- list(
	pvalues.out.noconsec,
	pvalues.cov.out.noconsec,
	pvalues.match.out.noconsec,
	pvalues.match.cov.out.noconsec,
	pvalues.trends.out.noconsec,
	pvalues.trendsSqrd.out.noconsec
	)

successes <- c(
	sum(diff_dat_no_consec$cp),
	Nsucc.cov.out.noconsec,
	Nsucc.match.out.noconsec,
	Nsucc.match.cov.out.noconsec,
	Nsucc.cov.trend.out.noconsec,
	Nsucc.trendssqr.out.noconsec
	)

attempts <- c(
	nrow(diff_dat_no_consec),
	N.cov.out.noconsec,
	N.match.out.noconsec,
	N.match.cov.out.noconsec,
	N.trend.out.noconsec,
	N.trendssqr.out.noconsec
	)

TableA2 <- stargazer(
	did.out.noconsec,
	estims.cov.noconsec,
	estims.match.noconsec,
	estims.match.cov.noconsec,
	did.out_trend.noconsec,
	estims.trendsqrd.noconsec,
	coef = coefsList,
	se = seList,
	p = pList,
	p.auto = FALSE,
	omit = c('factor', 'region', 'Constant'),
	dep.var.caption  = "\\emph{Dependent variable}: First Difference in State Repression",
	dep.var.label  = c('', '', '', '', '', ''),
	column.labels   = c('Baseline Models', "With CBPS Weighing", 'With Trends'),
	column.separate = c(2, 2, 2),
	dep.var.labels.include  = FALSE,
	align = TRUE,
	order = c('^cp', '^cheibubMil', '^ythblgap', '^cheibubDemo', '^ln', '^CW', '^counterbalancing', '^ncp', '^ncpFl', '^trend', '^trends'),
	omit.stat = c('rsq', 'f', 'ser', 'adj.rsq', 'n'),
	add.lines = list(
		c("CBPS Weighting", '$No$', '$No$', '$Yes$', '$Yes$', '$No$', '$No$'),#, '$No$', '$No$', '$Yes$', '$Yes$'),
		c("Observations", attempts),#, '$No$', '$No$', '$Yes$', '$Yes$'),
		c("Successful Coups", successes)#, '$No$', '$No$', '$Yes$', '$Yes$'),
		),
	notes.append = FALSE,
	float = FALSE
	)

#--------------------------------------
#Table A3
#--------------------------------------

did_ctrends_out_dat <- na.omit(diff_dat[, c('DV_diff', 'cp', 'nAttempts', 'ccode', 'Year', 'trend.pt1', 'trendsqrd.pt1', 'lm_weights')])
did_ctrends_out_dat$ccode <- factor(as.character(did_ctrends_out_dat$ccode))
did_ctrends_out_dat$ID <- 1:nrow(did_ctrends_out_dat)

did.ctrends.out.boot <- boot_lm(dat = did_ctrends_out_dat, DV_diff ~ cp + factor(nAttempts) + trend.pt1*factor(ccode), cluster_var = 'ID', strata = 'nAttempts')
did.ctrends.out <- lm(DV_diff ~ cp + factor(nAttempts) + trend.pt1*factor(ccode), data = did_ctrends_out_dat, y = TRUE, weights = lm_weights)

coefs.ctrends.out <- coef(did.ctrends.out)
SEs.ctrends.out <- apply(did.ctrends.out.boot, 2, function(x) sd(x, na.rm=T))
pvalues.ctrends.out <- 2*pnorm(-abs(coefs.ctrends.out/SEs.ctrends.out))

#--------------------
did.ctrendsSqrd.out.dat <- na.omit(diff_dat[, c('DV_diff', 'cp', 'nAttempts', 'ccode', 'Year', 'trend.pt1', 'trendsqrd.pt1', 'lm_weights')])
did.ctrendsSqrd.out.dat$ccode <- factor(as.character(did.ctrendsSqrd.out.dat$ccode))
did.ctrendsSqrd.out.dat$ID <- 1:nrow(did.ctrendsSqrd.out.dat)

did.ctrends.out.boot <- boot_lm(dat = did_ctrends_out_dat, DV_diff ~ cp + factor(nAttempts) + trend.pt1*factor(ccode) + trendsqrd.pt1*factor(ccode), cluster_var = 'ID', strata = 'nAttempts')
did.ctrendsSqrd.out <- lm(DV_diff ~ cp + factor(nAttempts) + trend.pt1*factor(ccode) + trendsqrd.pt1*factor(ccode), data = did.ctrendsSqrd.out.dat, y = TRUE, weights = lm_weights)

coefs.ctrendsSqrd.out <- coef(did.ctrends.out)
SEs.ctrendsSqrd.out <- apply(did.ctrends.out.boot, 2, function(x) sd(x, na.rm=T))
pvalues.ctrendsSqrd.out <- 2*pnorm(-abs(coefs.ctrendsSqrd.out/SEs.ctrendsSqrd.out))

pList <- list(
	pvalues.ctrends.out,
	pvalues.ctrendsSqrd.out
	)
seList <- list(
	SEs.ctrends.out,
	SEs.ctrendsSqrd.out
	)
successes <- c(
	sum(did_ctrends_out_dat$cp),
	sum(did.ctrendsSqrd.out.dat$cp)
	)
attempts <- c(
	nrow(did_ctrends_out_dat),
	nrow(did.ctrendsSqrd.out.dat)
	)

TableA3 <- stargazer(
	did.ctrends.out,#,
	did.ctrendsSqrd.out,
	se = seList,
	p = pList,
	p.auto = FALSE,
	omit = c('factor', 'region', 'trend', 'Constant'),
	dep.var.caption  = "\\emph{Dependent variable}: First Difference in State Repression",
	dep.var.label  = c('', ''),
	column.labels   = c('Country Trends', 'Country Trends Squared'),
	dep.var.labels.include  = FALSE,
	align = TRUE,
	omit.stat = c('rsq', 'f', 'ser', 'adj.rsq', 'n'),
	add.lines = list(
		c("Observations", attempts),
		c("Successful Coups", successes),
		c("Country-Specific Trends", 'Yes', 'Yes'),
		c("Country-Specific Trends Squared", 'No', 'Yes')
		),
	notes.append = FALSE,
	float = FALSE
	)

#--------------------------------------
#Figure A1
#--------------------------------------
D_bsc <- read.csv('ts_bsc.csv')
bsc.keep <- !is.na(D_bsc$cpAtL1)&!is.na(D_bsc$latentmean)&!is.na(D_bsc$latentmeanL2)
dat.bsc_t <- D_bsc[bsc.keep,]
formul.B <- latentmean ~ cpAtL1*latent.pt1 + factor(Year) + trend*factor(ccode) + trendsqrd*factor(ccode) + trendcube*factor(ccode) # + trendfrth*factor(ccode) + shrYth1524
dat.bsc_t$ccode <- as.factor(as.character(dat.bsc_t$ccode))

coefs <- NULL
ses <- NULL
for (i in 1:10) {
	dat.bsc_t$latent.pt1 <- dat.bsc_t[[paste('latentmean_post_', i, '_L2', sep = '')]]
	dat.bsc_t <- dat.bsc_t[!is.na(dat.bsc_t$latent.pt1), ]
	dat.bsc_t$ccode <- factor(dat.bsc_t$ccode)
	heter.out <- panel_reg(dat.bsc_t, update(formul.B,  . ~ . + latent.pt1*cpAtL1), cus_weights = NULL )
	coefs <- rbind(coefs, heter.out$coefs[grepl('(^cp|^latent)', rownames(heter.out$coefs)),][,1])
	ses <- rbind(ses, heter.out$coefs[grepl('(^cp|^latent)', rownames(heter.out$coefs)),][,2])
	print(i)
}
c_pth <- as.vector(mi.meld(coefs, ses)$q.mi)
s_pth <- as.vector(mi.meld(coefs, ses)$se.mi)
vrNms <- colnames(mi.meld(coefs, ses)$q.mi)
names(c_pth) <- vrNms
names(s_pth) <- vrNms

rnms1 <- which(vrNms == 'cpAtL1')
rnms2 <- which(vrNms == 'cpAtL1:latent.pt1')
pld_cp_est <- c_pth[c(rnms1, rnms2)]
pld_cp_var <- s_pth[c(rnms1, rnms2)]
pool_varVcov <- matrix(nrow = 2, ncol = 2, 0)
diag(pool_varVcov) <- (pld_cp_var)^2
draws <- mvrnorm(n = 100000, mu = pld_cp_est, Sigma = pool_varVcov)
colnames(draws) <- c('cpAtL1', 'cpAtL1:latent.pt1')

latentmean_levels <-  seq(-4, 4, 0.01)
frstTerm <- as.matrix(data.frame(rep(1, length(latentmean_levels)), latentmean_levels))
marginal_effects <- frstTerm %*% t(draws[, grepl('^cp', colnames(draws))])
means <- apply(marginal_effects, 1, mean)
l90 <- apply(marginal_effects, 1, function(x) quantile(x, 0.05))
u90 <- apply(marginal_effects, 1, function(x) quantile(x, 0.95))
l95 <- apply(marginal_effects, 1, function(x) quantile(x, 0.025))
u95 <- apply(marginal_effects, 1, function(x) quantile(x, 0.975))
l99 <- apply(marginal_effects, 1, function(x) quantile(x, 0.005))
u99 <- apply(marginal_effects, 1, function(x) quantile(x, 0.995))

pdf('FigureA1.pdf', width = 16*0.5, height = 12*0.5)
	attach(mtcars)
	plot(x = latentmean_levels,
		y = means,
		type = 'l',
		xlim = c(-2, 3),
		ylim = c(-0.15, 0.33),
		xlab = 'Pre-Treatment Level of State Repression',
		ylab = 'Marginal Effect of Coup Attempt on State Repression',
		xaxt = 'n',
		yaxt = 'n',
		lty = 1,
		lwd = 1
		)
	axis(1, at = -2:3, labels = -2:3)
	axis(2, at = seq(-0.2, 0.3, by = 0.1), labels = seq(-0.2, 0.3, by = 0.1), las = 1)
	polygon(c(latentmean_levels, rev(latentmean_levels)),
		c((l90), rev((u90))),
		border = FALSE,
		fillOddEven = TRUE,
		col = rgb(171/255, 171/255, 171/255, alpha= 0.55))
	polygon(c(latentmean_levels, rev(latentmean_levels)),
		c((l95), rev((u95))),
		border = FALSE,
		fillOddEven = TRUE,
		col = rgb(171/255, 171/255, 171/255, alpha= 0.25))
	abline(0, 0, lty =1, lwd = 0.5)
	abline(v = 0, lty = 2, col = 'grey')
	abline(v = 1, lty = 2, col = 'grey')
	abline(v = 2, lty = 2, col = 'grey')
	abline(v = 3, lty = 2, col = 'grey')
	abline(v = -1, lty = 2, col = 'grey')
	abline(v = -2, lty = 2, col = 'grey')
	abline(v = -3, lty = 2, col = 'grey')
dev.off()


#---------------------------------------------
#Figure A2
#---------------------------------------------
D_pth <- read.csv('ts_pth.csv')
formul <- latentmean ~ cpLead2 + cpLead1 + cp + cpL1 + cpL2 + cpL3 + cpL4 + cpL5 + factor(ccode) + trend*factor(ccode) +factor(Year)+ trendsqrd*factor(ccode) + trendcube*factor(ccode) # + cheibubDemoL2
minYear <- 1950; maxYear <- 2013 #######                  #####

placeb.dat.pth <- na.omit(D_pth[D_pth$Year >= minYear & D_pth$Year <= maxYear, c('ccode', all.vars(formul))])
placeb.dat.rch <- na.omit(D_rch[D_rch$Year >= minYear & D_rch$Year <= maxYear, c('ccode', all.vars(formul))])

placeb.dat.pth$ccode <- factor(as.character(placeb.dat.pth$ccode))
placeb.dat.rch$ccode <- factor(as.character(placeb.dat.rch$ccode))

placeb.reg.pth <- lm(formul, data = placeb.dat.pth)
placeb.reg.rch <- lm(formul, data = placeb.dat.rch)

CovVar.pth <- vcovCluster(model = placeb.reg.pth, cluster = placeb.dat.pth$ccode)
CovVar.rch <- vcovCluster(model = placeb.reg.rch, cluster = placeb.dat.rch$ccode)

results.pth <- coeftest(placeb.reg.pth, CovVar.pth, df.residual(placeb.reg.rch))
results.rch <- coeftest(placeb.reg.rch, CovVar.rch, df.residual(placeb.reg.rch))

Means.pth <- c(
	results.pth[,1][names(results.pth[,1]) == 'cpLead2'],
	results.pth[,1][names(results.pth[,1]) == 'cpLead1'],
	results.pth[,1][names(results.pth[,1]) == 'cp'],
	results.pth[,1][names(results.pth[,1]) == 'cpL1'],
	results.pth[,1][names(results.pth[,1]) == 'cpL2'],
	results.pth[,1][names(results.pth[,1]) == 'cpL3'],
	results.pth[,1][names(results.pth[,1]) == 'cpL4'],
	results.pth[,1][names(results.pth[,1]) == 'cpL5']
	)
Means.rch <- c(
	results.rch[,1][names(results.rch[,1]) == 'cpLead2'],
	results.rch[,1][names(results.rch[,1]) == 'cpLead1'],
	results.rch[,1][names(results.rch[,1]) == 'cp'],
	results.rch[,1][names(results.rch[,1]) == 'cpL1'],
	results.rch[,1][names(results.rch[,1]) == 'cpL2'],
	results.rch[,1][names(results.rch[,1]) == 'cpL3'],
	results.rch[,1][names(results.rch[,1]) == 'cpL4'],
	results.rch[,1][names(results.rch[,1]) == 'cpL5']
	)

SEs.pth <- c(
	results.pth[,2][names(results.pth[,2]) == 'cpLead2'],
	results.pth[,2][names(results.pth[,2]) == 'cpLead1'],
	results.pth[,2][names(results.pth[,2]) == 'cp'],
	results.pth[,2][names(results.pth[,2]) == 'cpL1'],
	results.pth[,2][names(results.pth[,2]) == 'cpL2'],
	results.pth[,2][names(results.pth[,2]) == 'cpL3'],
	results.pth[,2][names(results.pth[,2]) == 'cpL4'],
	results.pth[,2][names(results.pth[,2]) == 'cpL5']
	)

SEs.rch <- c(
	results.rch[,2][names(results.rch[,2]) == 'cpLead2'],
	results.rch[,2][names(results.rch[,2]) == 'cpLead1'],
	results.rch[,2][names(results.rch[,2]) == 'cp'],
	results.rch[,2][names(results.rch[,2]) == 'cpL1'],
	results.rch[,2][names(results.rch[,2]) == 'cpL2'],
	results.rch[,2][names(results.rch[,2]) == 'cpL3'],
	results.rch[,2][names(results.rch[,2]) == 'cpL4'],
	results.rch[,2][names(results.rch[,2]) == 'cpL5']
	)

l.pth <- qnorm(p = 0.005, mean = Means.pth, sd = SEs.pth)
u.pth <- qnorm(p = 0.995, mean = Means.pth, sd = SEs.pth)

l.rch <- qnorm(p = 0.005, mean = Means.rch, sd = SEs.rch)
u.rch <- qnorm(p = 0.995, mean = Means.rch, sd = SEs.rch)

ylim = c(-0.05, 0.31)
pdf('FigureA2.pdf', width = 11, height = 8)
	par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(1.9,3.6,1.1,1.1))
	par(mfrow=c(2,1))
	plot(x = -2:5, y = Means.pth, xlim = c(-2.5, 5.5), ylim = ylim, xaxt = "n", yaxt = "n", ylab = '', las = 2, xlab = '')
	lines(x = -2:5, y = Means.pth, lty = 3)
	abline(h = 0, lty = 2)
	abline(v = 0, lty = 2)
	arrows(x0 = -2:5, y0 = l.pth, y1 = u.pth, angle = 90, code = 0, lwd=2)
	title('a) Powell and Thyne Coup Data')
	axis(side = 1, at = -2:5, labels = c('2 Years Prior', '1 Year Prior', 'Year of Event', '1 Year After', '2 Year After', '3 Year After', '4 Years After', '5 Years After'))
	title(ylab="Effect on State Repression", line=4, cex.lab=1.2)
	axis(side = 1, at = -2:5, labels = c('2 Years Prior', '1 Year Prior', 'Year of Event', '1 Year After', '2 Year After', '3 Year After', '4 Years After', '5 Years After'), cex.lab=4.2)
	axis(side = 2, at = seq(-0.1, 0.5, by = .1), labels = seq(-0.1, 0.5, by = .1), cex.lab=2.2, las = 2)

	plot(x = -2:5, y = Means.rch, xlim = c(-2.5, 5.5), ylim = ylim, xaxt = "n", yaxt = "n", ylab = '', xlab = '')#, ylab = 'Estimated Effect on State Repression', las = 2)#, xlab = 'Time Passage Relative to Coup Attempt')
	lines(x = -2:5, y = Means.rch, lty = 3)
	abline(h = 0, lty = 2)
	abline(v = 0, lty = 2)
	arrows(x0 = -2:5, y0 = l.rch, y1 = u.rch, angle = 90, code = 0, lwd=2)
	title('b) Archigos Coup Data')
	axis(side = 1, at = -2:5, labels = c('2 Years Prior', '1 Year Prior', 'Year of Event', '1 Year After', '2 Year After', '3 Year After', '4 Years After', '5 Years After'))
	title(xlab="Time Passage Relative to Successful Coup", line=3, cex.lab=1.2)
	title(ylab="Effect on State Repression", line=4, cex.lab=1.2)
	axis(side = 1, at = -2:5, labels = c('2 Years Prior', '1 Year Prior', 'Year of Event', '1 Year After', '2 Year After', '3 Year After', '4 Years After', '5 Years After'), cex.lab=4.2)
	axis(side = 2, at = seq(-0.1, 0.5, by = .1), labels = seq(-0.1, 0.5, by = .1), cex.lab=2.2, las = 2)

dev.off()



#------------------------------
#Figure A3
#------------------------------

formul <- DV_diff ~ cp + factor(nAttempts)
dat.cheibmil <- na.omit(diff_dat[diff_dat$cheibubMil.pt1 %in% 1 & diff_dat$cheibubDemo.pt1 %in% 0, c('DV_diff', 'cp', 'nAttempts', 'ccode', 'Year')])
dat.cheibmil$ID <- 1:nrow(dat.cheibmil)
dat.cheibmil$nAttempts[dat.cheibmil$nAttempts>=2] <- 2
dat.cheibmil$nAttempts <-  factor(as.character(dat.cheibmil$nAttempts))

diff_temp1 <- subset(dat.cheibmil, nAttempts ==1)
diff_temp2 <- subset(dat.cheibmil, nAttempts ==2)
pscore1 <- mean(diff_temp1$cp)
pscore2 <- mean(diff_temp2$cp)
dat.cheibmil$lm_weights[dat.cheibmil$nAttempts == 1 & dat.cheibmil$cp == 1] <- 1/pscore1
dat.cheibmil$lm_weights[dat.cheibmil$nAttempts == 1 & dat.cheibmil$cp == 0] <- 1/(1-pscore1)
dat.cheibmil$lm_weights[dat.cheibmil$nAttempts == 2 & dat.cheibmil$cp == 1] <- 1/pscore2
dat.cheibmil$lm_weights[dat.cheibmil$nAttempts == 2 & dat.cheibmil$cp == 0] <- 1/(1-pscore2)

cheibmil_boot_AR <- boot_lm(dat = dat.cheibmil, formul = formul, cluster_var = 'ID', strata = 'nAttempts')
l90_mil <- quantile(cheibmil_boot_AR[,colnames(cheibmil_boot_AR)=='cp'], 0.05)
u90_mil <- quantile(cheibmil_boot_AR[,colnames(cheibmil_boot_AR)=='cp'], 0.95)
did.out.cheibmil <- lm(DV_diff ~ cp + factor(nAttempts),  dat.cheibmil, weights = lm_weights)
coef_mil <- coef(did.out.cheibmil)[names(coef(did.out.cheibmil)) == 'cp']

#Non-Mil Regimes
dat.cheib.non.mil <- na.omit(diff_dat[diff_dat$cheibubMil.pt1 %in% 0 & diff_dat$cheibubDemo.pt1 %in% 0, c('DV_diff', 'cp', 'nAttempts', 'ccode', 'Year')])
dat.cheib.non.mil$ID <- 1:nrow(dat.cheib.non.mil)
dat.cheib.non.mil$nAttempts[dat.cheib.non.mil$nAttempts>=2] <- 2
dat.cheib.non.mil$nAttempts <-  factor(as.character(dat.cheib.non.mil$nAttempts))

diff_temp1 <- subset(dat.cheib.non.mil, nAttempts ==1)
diff_temp2 <- subset(dat.cheib.non.mil, nAttempts ==2)
pscore1 <- mean(diff_temp1$cp)
pscore2 <- mean(diff_temp2$cp)
dat.cheib.non.mil$lm_weights[dat.cheib.non.mil$nAttempts == 1 & dat.cheib.non.mil$cp == 1] <- 1/pscore1
dat.cheib.non.mil$lm_weights[dat.cheib.non.mil$nAttempts == 1 & dat.cheib.non.mil$cp == 0] <- 1/(1-pscore1)
dat.cheib.non.mil$lm_weights[dat.cheib.non.mil$nAttempts == 2 & dat.cheib.non.mil$cp == 1] <- 1/pscore2
dat.cheib.non.mil$lm_weights[dat.cheib.non.mil$nAttempts == 2 & dat.cheib.non.mil$cp == 0] <- 1/(1-pscore2)

cheinonbmil_boot_AR <- boot_lm(dat = dat.cheib.non.mil, formul = formul, cluster_var = 'ID', strata = 'nAttempts')
l90_nonmil <- quantile(cheinonbmil_boot_AR[,colnames(cheinonbmil_boot_AR)=='cp'], 0.05)
u90_nonmil <- quantile(cheinonbmil_boot_AR[,colnames(cheinonbmil_boot_AR)=='cp'], 0.95)
did.out.cheibnonmil <- lm(DV_diff ~ cp + factor(nAttempts),  dat.cheib.non.mil, weights = lm_weights)
coef_nonmil <- coef(did.out.cheibnonmil)[names(coef(did.out.cheibnonmil)) == 'cp']

l90 <- c(
	l90_mil,
	l90_nonmil
	)
u90 <- c(
	u90_mil,
	u90_nonmil
	)
coefs <- c(
	coef_mil,
	coef_nonmil
	)

xlabels1 <- c(
	'Military Regime',
	'Non-Military Regime'
	)
xlabels2 <- c(
	paste('N = ', nrow(dat.cheibmil)),
	paste('N = ', nrow(dat.cheib.non.mil))
	)

pdf('FigureA3.pdf', width = 6, height = 4)
	par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.9,3.6,1.1,1.1))
	plot(x = 1:2, y = coefs, xlim = c(0.5, 2.5), ylim = c(-0.1, 0.45), xaxt = "n", yaxt = "n", ylab = 'Treatment Effect', las = 2, xlab = '')
	axis(side = 1, at = c(1, 2), labels = xlabels1)
	axis(side = 1, at = c(1, 2), labels = xlabels2, line = 1.3, tick = FALSE)
	axis(side = 2, at = c(0.4, 0.3, 0.2, 0.1, 0, -0.1, -0.2), labels = c(0.4, 0.3, 0.2, 0.1, 0, -0.1, -0.2), las = 3)
	abline(h = 0, lty = 2, lwd = 0.5)
	arrows(x0 = 1:2, y0 = l90, y1 = u90, angle = 90, code = 3)
dev.off()


#------------------------------
#Figure A4
#------------------------------

dat.mil <- na.omit(diff_dat[diff_dat$gwf_military.pt1 %in% 1, c('DV_diff', 'cp', 'nAttempts', 'ccode', 'Year')])
dat.mil$ID <- 1:nrow(dat.mil)
dat.mil$nAttempts[dat.mil$nAttempts>=2] <- 2
dat.mil$nAttempts <- factor(as.character(dat.mil$nAttempts))

diff_temp1 <- subset(dat.mil, nAttempts ==1)
diff_temp2 <- subset(dat.mil, nAttempts ==2)
pscore1 <- mean(diff_temp1$cp)
pscore2 <- mean(diff_temp2$cp)
dat.mil$lm_weights[dat.mil$nAttempts == 1 & dat.mil$cp == 1] <- 1/pscore1
dat.mil$lm_weights[dat.mil$nAttempts == 1 & dat.mil$cp == 0] <- 1/(1-pscore1)
dat.mil$lm_weights[dat.mil$nAttempts == 2 & dat.mil$cp == 1] <- 1/pscore2
dat.mil$lm_weights[dat.mil$nAttempts == 2 & dat.mil$cp == 0] <- 1/(1-pscore2)

mil_boot_AR <- boot_lm(dat = dat.mil, formul = formul, cluster_var = 'ID', strata = 'nAttempts')
l90_mil <- quantile(mil_boot_AR[,colnames(mil_boot_AR)=='cp'], 0.05)
u90_mil <- quantile(mil_boot_AR[,colnames(mil_boot_AR)=='cp'], 0.95)
did.out.mil <- lm(DV_diff ~ cp + factor(nAttempts),  dat.mil, weights = lm_weights)
coef_mil <- coef(did.out.mil)[names(coef(did.out.mil)) == 'cp']
N.mil <- nrow(dat.mil)


#GWF Party
dat.part <- na.omit(diff_dat[diff_dat$gwf_party.pt1 %in% 1, c('DV_diff', 'cp', 'nAttempts', 'ccode', 'Year')])
dat.part$ID <- 1:nrow(dat.part)
dat.part$nAttempts[dat.part$nAttempts>=2] <- 2
dat.part$nAttempts <- factor(as.character(dat.part$nAttempts))

diff_temp1 <- subset(dat.part, nAttempts ==1)
diff_temp2 <- subset(dat.part, nAttempts ==2)
pscore1 <- mean(diff_temp1$cp)
pscore2 <- mean(diff_temp2$cp)
dat.part$lm_weights[dat.part$nAttempts == 1 & dat.part$cp == 1] <- 1/pscore1
dat.part$lm_weights[dat.part$nAttempts == 1 & dat.part$cp == 0] <- 1/(1-pscore1)
dat.part$lm_weights[dat.part$nAttempts == 2 & dat.part$cp == 1] <- 1/pscore2
dat.part$lm_weights[dat.part$nAttempts == 2 & dat.part$cp == 0] <- 1/(1-pscore2)

part_boot_AR <- boot_lm(dat = dat.part, formul = formul, cluster_var = 'ID', strata = 'nAttempts')
l90_part <- quantile(part_boot_AR[,colnames(part_boot_AR)=='cp'], 0.05)
l90_part
u90_part <- quantile(part_boot_AR[,colnames(part_boot_AR)=='cp'], 0.95)
u90_part
did.out.part <- lm(DV_diff ~ cp + factor(nAttempts),  dat.part, weights = lm_weights)
coef_part <- coef(did.out.part)[names(coef(did.out.part)) == 'cp']
N.part <- nrow(dat.part)

#GWF Pers
dat.pers <- na.omit(diff_dat[diff_dat$gwf_personal.pt1 %in% 1, c('DV_diff', 'cp', 'nAttempts', 'ccode', 'Year')])
dat.pers$ID <- 1:nrow(dat.pers)
dat.pers$nAttempts[dat.pers$nAttempts>=2] <- 2
dat.pers$nAttempts <- factor(dat.pers$nAttempts)

diff_temp1 <- subset(dat.pers, nAttempts ==1)
diff_temp2 <- subset(dat.pers, nAttempts ==2)
pscore1 <- mean(diff_temp1$cp)
pscore2 <- mean(diff_temp2$cp)
dat.pers$lm_weights[dat.pers$nAttempts == 1 & dat.pers$cp == 1] <- 1/pscore1
dat.pers$lm_weights[dat.pers$nAttempts == 1 & dat.pers$cp == 0] <- 1/(1-pscore1)
dat.pers$lm_weights[dat.pers$nAttempts == 2 & dat.pers$cp == 1] <- 1/pscore2
dat.pers$lm_weights[dat.pers$nAttempts == 2 & dat.pers$cp == 0] <- 1/(1-pscore2)

pers_boot_AR <- boot_lm(dat = dat.pers, formul = formul, cluster_var = 'ID', strata = 'nAttempts')
l90_pers <- quantile(pers_boot_AR[,colnames(pers_boot_AR)=='cp'], 0.05)
u90_pers <- quantile(pers_boot_AR[,colnames(pers_boot_AR)=='cp'], 0.95)
did.out.path <- lm(DV_diff ~ cp + factor(nAttempts),  dat.pers, weights = lm_weights)
coef_pers <- coef(did.out.path)[names(coef(did.out.path)) == 'cp']
N.pers <- nrow(dat.pers)

#GWF Monarchy
dat.mon <- na.omit(diff_dat[diff_dat$gwf_monarch.pt1 %in% 1, c('DV_diff', 'cp', 'nAttempts', 'ccode', 'Year')])
dat.mon$ID <- 1:nrow(dat.mon)
dat.mon$nAttempts[dat.mon$nAttempts>=2] <- 2
diff_temp1 <- subset(dat.mon, nAttempts ==1)
diff_temp2 <- subset(dat.mon, nAttempts ==2)
pscore1 <- mean(diff_temp1$cp)
pscore2 <- mean(diff_temp2$cp)
dat.mon$lm_weights[dat.mon$nAttempts == 1 & dat.mon$cp == 1] <- 1/pscore1
dat.mon$lm_weights[dat.mon$nAttempts == 1 & dat.mon$cp == 0] <- 1/(1-pscore1)
dat.mon$lm_weights[dat.mon$nAttempts == 2 & dat.mon$cp == 1] <- 1/pscore2
dat.mon$lm_weights[dat.mon$nAttempts == 2 & dat.mon$cp == 0] <- 1/(1-pscore2)

did.out.mon <- lm(DV_diff ~ cp + factor(nAttempts),  dat.mon, weights = lm_weights)
l90_mon <- coef(did.out.mon)[2] - 1.645*summary(did.out.mon)$coefficients[2 ,2]
u90_mon <- coef(did.out.mon)[2] + 1.645*summary(did.out.mon)$coefficients[2 ,2]
coef_mon <- coef(did.out.mon)[names(coef(did.out.mon)) == 'cp']
N.mon <- nrow(dat.mon)

l90 <- c(l90_mil, l90_part, l90_pers, l90_mon)
u90 <- c(u90_mil, u90_part, u90_pers, u90_mon)
coefs <- c(coef_mil, coef_part, coef_pers, coef_mon)

xlabels1 <- c(
	'Mil. Reg.',
	'Part. Reg.',
	'Pers. Reg.',
	'Monarch.'
	)
xlabels2 <- c(
	paste('N = ', N.mil),
	paste('N = ', N.part),
	paste('N = ', N.pers),
	paste('N = ', N.mon)
	)

pdf('FigureA4.pdf', width = 6, height = 4)
	par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.9,3.6,1.1,1.1))
	plot(x = 1:4, y = coefs, xlim = c(0.5, 4.5), ylim = c(-.3, .8), xaxt = "n", yaxt = "n", ylab = 'Treatment Effect', las = 2, xlab = '')
	axis(side = 1, at = c(1, 2, 3, 4), labels = xlabels1)
	axis(side = 1, at = c(1, 2, 3, 4), labels = xlabels2, line = 1.3, tick = FALSE)
	axis(side = 2, at = c(0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0, -0.1, -0.2), labels = c(0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0, -0.1, -0.2), las = 3)
	abline(h = 0, lty = 2, lwd = 0.5)
	arrows(x0 = 1:4, y0 = l90, y1 = u90, angle = 90, code = 3)
dev.off()

#------------------------------
#Section B1
#------------------------------
library(cobalt)
#Calculation of standardized means difference for Powell and Thyne, differenced dataset
balanceVars.pth <- balanceVars[!(balanceVars %in% 'ncp.bsc')]
balanceVars.pth <- balanceVars.pth[!(balanceVars.pth %in% 'ncpFl.bsc')]
balanceVars.pth <- c(balanceVars.pth, 'ncp', 'ncpFl')
balanceVars.pt1 <- gsub('L1', '.pt1', balanceVars.pth)
balanceVars.pt1 <- gsub('^ncpFl$', 'ncpFl.pt1', balanceVars.pt1)
balanceVars.pt1 <- gsub('^ncp$', 'ncp.pt1', balanceVars.pt1)
balanceVars.pt1 <- gsub('CW', 'CW.pt1', balanceVars.pt1)

identify.binary.vars <- apply(diff_dat_pth[, balanceVars.pt1], 2, checkBinaryTrait)
sddMnsdff_pth <- NULL
balanceTable <- matrix("", nrow = 0 * length(balanceVars.pt1), ncol = 5)
for (var in balanceVars.pt1){
	print (var)
	dat_pth_t <- na.omit(diff_dat_pth[, c('cp', var)])
	frml <- as.formula(paste('cp~', var))
	a <- CBPS(frml, data = dat_pth_t, m.threshold = 0.25)
	sddMnsdff_pth_t <- bal.tab(a, un = TRUE)$Balance$Diff.Un[2]
	sddMnsdff_pth <- c(sddMnsdff_pth, sddMnsdff_pth_t)
}

names(sddMnsdff_pth) <- balanceVars.pt1
abs(sddMnsdff_pth)[order(abs(sddMnsdff_pth))]


#------------------------------
#Table B1
#------------------------------
controls.pth <- c(
				'cheibubMil',
				'ythblgap',
				'cheibubDemo',
				'ln.rCOWmilexpc',
				'CW',
				'counterbalancing',
				'ncp',
				'ncpFl'
			)
match.on.pth <- c(
	 'e_migdppcln',
	 'cheibubLegis',
	 'ln.milperpc',
	 'counterbalancing',
	 'oil_binary',
	 'purges',
	 'ncp',
	 'ncpFl',
	 'OppNonViolBanks',
	 'ethfrac'
	 )

diff_dat_pth$nAttempts[as.integer(as.character(diff_dat_pth$nAttempts)) >= 3] <- 3
diff_temp1 <- subset(diff_dat_pth, nAttempts ==1)
diff_temp2 <- subset(diff_dat_pth, nAttempts ==2)
diff_temp3 <- subset(diff_dat_pth, nAttempts ==3)
pscore1 <- mean(diff_temp1$cp)
pscore2 <- mean(diff_temp2$cp)
pscore3 <- mean(diff_temp3$cp)

diff_dat_pth$lm_weights[as.integer(as.character(diff_dat_pth$nAttempts)) == 1 & diff_dat_pth$cp == 1] <- 1/pscore1
diff_dat_pth$lm_weights[as.integer(as.character(diff_dat_pth$nAttempts)) == 1 & diff_dat_pth$cp == 0] <- 1/(1-pscore1)
diff_dat_pth$lm_weights[as.integer(as.character(diff_dat_pth$nAttempts)) == 2 & diff_dat_pth$cp == 1] <- 1/pscore2
diff_dat_pth$lm_weights[as.integer(as.character(diff_dat_pth$nAttempts)) == 2 & diff_dat_pth$cp == 0] <- 1/(1-pscore2)
diff_dat_pth$lm_weights[as.integer(as.character(diff_dat_pth$nAttempts)) == 3 & diff_dat_pth$cp == 1] <- 1/pscore3
diff_dat_pth$lm_weights[as.integer(as.character(diff_dat_pth$nAttempts)) == 3 & diff_dat_pth$cp == 0] <- 1/(1-pscore3)

diff_dat_pth$ID <- 1:nrow(diff_dat_pth)

diff_dat_pth_mi$ID <- 1:nrow(diff_dat_pth_mi)
diff_dat_pth_mi$ccode <- factor(as.character(diff_dat_pth_mi$ccode))
diff_dat_pth_mi$w <- NA
match.formul <- formula(paste('cp ~ ', paste(paste(match.on.pth, '.pt1', sep = '', collapse = ' + '), sep = ''), ' + nAttempts', sep = ''))
for (i in levels(diff_dat_pth_mi$.imp)) {
	in.imp <- diff_dat_pth_mi$.imp == i
	m.out <- CBPS(match.formul, standardize = F, ATT = 0, data = diff_dat_pth_mi[in.imp,])
	diff_dat_pth_mi$w[in.imp] <- m.out$weights
	diff_dat_pth_mi$ID[in.imp] <- diff_dat_pth$ID
	diff_dat_pth_mi$lm_weights[in.imp] <- diff_dat_pth$lm_weights
}

#Baseline
#------------------------------------------
library(lmtest)
formul <- DV_diff ~ cp + factor(nAttempts)
did.out_boot.pth <- boot_lm(dat = diff_dat_pth, formul = formul, cluster_var = 'ID', strata = 'nAttempts')
did.out.pth <- lm(formul, data = diff_dat_pth, weights = lm_weights)
coefs.out.pth <- coef(did.out.pth)
SEs.out.pth <- apply(did.out_boot.pth, 2, sd)
pvalues.out.pth <- 2*pnorm(-abs(coefs.out.pth/SEs.out.pth))

estims.cov.pth <- lm_mi(dat_t = diff_dat_pth_mi, controls = controls.pth, reweigh = 'lm_weights')
coefs.cov.out.pth <- estims.cov.pth$coefs[,1]
SEs.cov.out.pth <- estims.cov.pth$coefs[,2]
pvalues.cov.out.pth <- estims.cov.pth$coefs[,3]
N.cov.out.pth <- estims.cov.pth$N
Nsucc.cov.out.pth <- estims.cov.pth$Nsucc

#Matching
#------------------------------------------
estims.match.pth <- lm_mi(dat_t = diff_dat_pth_mi, controls = NULL, reweigh = 'w')
coefs.match.out.pth <- estims.match.pth$coefs[,1]
SEs.match.out.pth <- estims.match.pth$coefs[,2]
pvalues.match.out.pth <- estims.match.pth$coefs[,3]
N.match.out.pth <- estims.match.pth$N
Nsucc.match.out.pth <- estims.match.pth$Nsucc

estims.match.cov.pth <- lm_mi(dat_t = diff_dat_pth_mi, controls = controls.pth, reweigh = 'w')
coefs.match.cov.out.pth <- estims.match.cov.pth$coefs[,1]
SEs.match.cov.out.pth <- estims.match.cov.pth$coefs[,2]
pvalues.match.cov.out.pth <- estims.match.cov.pth$coefs[,3]
N.match.cov.out.pth <- estims.match.cov.pth$N
Nsucc.match.cov.out.pth <- estims.match.cov.pth$Nsucc

#Trends
#------------------------------------------
formul_trend <- DV_diff ~ cp + factor(nAttempts) + trend.pt1

did.out_boot_trend.pth <- boot_lm(dat = diff_dat_pth, formul = formul_trend, cluster_var = 'ID', strata = 'nAttempts', R = 500)
did.out_trend.pth <- lm(formul_trend, data = diff_dat_pth, weights = lm_weights)
coefs.trends.out.pth <- coef(did.out_trend.pth)
SEs.trends.out.pth <- apply(did.out_boot_trend.pth, 2, sd)
pvalues.trends.out.pth <- 2*pnorm(-abs(coefs.trends.out.pth/SEs.trends.out.pth))
N.trend.out.pth <- nrow(diff_dat_pth)
Nsucc.cov.trend.out.pth <- sum(diff_dat_pth$cp)

estims.trendsqrd.pth <- lm_mi(dat_t = diff_dat_pth_mi, controls = c(controls.pth, 'trend', 'trendsqrd'), reweigh = 'lm_weights')
coefs.trendsSqrd.out.pth <- estims.trendsqrd.pth$coefs[,1]
SEs.trendsSqrd.out.pth <- estims.trendsqrd.pth$coefs[,2]
pvalues.trendsSqrd.out.pth <- estims.trendsqrd.pth$coefs[,3]
N.trendssqr.out.pth <- estims.trendsqrd.pth$N
Nsucc.trendssqr.out.pth <- estims.trendsqrd.pth$Nsucc

coefsList <- list(
	coefs.out.pth,
	coefs.cov.out.pth,
	coefs.match.out.pth,
	coefs.match.cov.out.pth,
	coefs.trends.out.pth,
	coefs.trendsSqrd.out.pth
	)

seList <- list(
	SEs.out.pth,
	SEs.cov.out.pth,
	SEs.match.out.pth,
	SEs.match.cov.out.pth,
	SEs.trends.out.pth,
	SEs.trendsSqrd.out.pth
	)

pList <- list(
	pvalues.out.pth,
	pvalues.cov.out.pth,
	pvalues.match.out.pth,
	pvalues.match.cov.out.pth,
	pvalues.trends.out.pth,
	pvalues.trendsSqrd.out.pth
	)

successes <- c(
	sum(diff_dat_pth$cp),
	Nsucc.cov.out.pth,
	Nsucc.match.out.pth,
	Nsucc.match.cov.out.pth,
	Nsucc.cov.trend.out.pth,
	Nsucc.trendssqr.out.pth
	)

attempts <- c(
	nrow(diff_dat_pth),
	N.cov.out.pth,
	N.match.out.pth,
	N.match.cov.out.pth,
	N.trend.out.pth,
	N.trendssqr.out.pth
	)

TableB1 <- stargazer(
	did.out.pth,
	estims.cov.pth$model,
	estims.match.pth$model,
	estims.match.cov.pth$model,
	did.out_trend.pth,
	estims.trendsqrd.pth$model,
	coef = coefsList,
	se = seList,
	p = pList,
	p.auto = FALSE,
	omit = c('factor', 'region', 'Constant'),
	dep.var.caption  = "\\emph{Dependent variable}: First Difference in State Repression",
	dep.var.label  = c('', '', '', '', '', ''),
	column.labels   = c('Baseline Models', "With CBPS Weighing", 'With Trends'),
	column.separate = c(2, 2, 2),
	dep.var.labels.include  = FALSE,
	align = TRUE,
	omit.stat = c('rsq', 'f', 'ser', 'adj.rsq', 'n'),
	add.lines = list(
		c("CBPS Weighing", '$No$', '$No$', '$Yes$', '$Yes$', '$No$', '$No$'),#, '$No$', '$No$', '$Yes$', '$Yes$'),
		c("Observations", attempts),#, '$No$', '$No$', '$Yes$', '$Yes$'),
		c("Successful Coups", successes)#, '$No$', '$No$', '$Yes$', '$Yes$'),
		),
	notes.append = FALSE,
	float = FALSE
	)

#------------------------------
#Figure B1
#------------------------------

dat_AR_pth <- subset(diff_dat_pth, cheibubDemo.pt1 == 0)
dat_AR_pth$ID <- 1:nrow(dat_AR_pth)

dat_AR_pth_n1 <- subset(dat_AR_pth, nAttempts == 1)
dat_AR_pth_n2 <- subset(dat_AR_pth, nAttempts == 2)
dat_AR_pth_n3 <- subset(dat_AR_pth, nAttempts == 3)
pscore1_AR <- mean(dat_AR_pth_n1$cp)
pscore2_AR <- mean(dat_AR_pth_n2$cp)
pscore3_AR <- mean(dat_AR_pth_n3$cp)

dat_AR_pth$lm_weights[dat_AR_pth$nAttempts == 1 & dat_AR_pth$cp == 1] <- 1/pscore1_AR
dat_AR_pth$lm_weights[dat_AR_pth$nAttempts == 1 & dat_AR_pth$cp == 0] <- 1/(1-pscore1_AR)
dat_AR_pth$lm_weights[dat_AR_pth$nAttempts == 2 & dat_AR_pth$cp == 1] <- 1/pscore2_AR
dat_AR_pth$lm_weights[dat_AR_pth$nAttempts == 2 & dat_AR_pth$cp == 0] <- 1/(1-pscore2_AR)
dat_AR_pth$lm_weights[dat_AR_pth$nAttempts == 3 & dat_AR_pth$cp == 1] <- 1/pscore3_AR
dat_AR_pth$lm_weights[dat_AR_pth$nAttempts == 3 & dat_AR_pth$cp == 0] <- 1/(1-pscore3_AR)

did.out_boot_AR_pth <- boot_lm(dat = dat_AR_pth, formul = DV_diff ~ cp + factor(nAttempts), cluster_var = 'ID', strata = 'nAttempts')
l90_autoc_pth <- quantile(did.out_boot_AR_pth[,colnames(did.out_boot_AR_pth)=='cp'], 0.05)
u90_autoc_pth <- quantile(did.out_boot_AR_pth[,colnames(did.out_boot_AR_pth)=='cp'], 0.95)

did.out.AR_pth <- lm(DV_diff ~ cp + factor(nAttempts),  dat_AR_pth, weights = lm_weights)
coef_AR_pth <- coef(did.out.AR_pth)[names(coef(did.out.AR_pth)) == 'cp']

dat_DEM_pth <- subset(diff_dat_pth, cheibubDemo.pt1 == 1)
dat_DEM_pth$nAttempts[dat_DEM_pth$nAttempts == 3] <- 2

dat_DEM_pth_n1 <- subset(dat_DEM_pth, nAttempts == 1)
dat_DEM_pth_n2 <- subset(dat_DEM_pth, nAttempts == 2)
pscore1_DEM <- mean(dat_DEM_pth_n1$cp)
pscore2_DEM <- mean(dat_DEM_pth_n2$cp)

dat_DEM_pth$lm_weights[dat_DEM_pth$nAttempts == 1 & dat_DEM_pth$cp == 1] <- 1/pscore1_DEM
dat_DEM_pth$lm_weights[dat_DEM_pth$nAttempts == 1 & dat_DEM_pth$cp == 0] <- 1/(1-pscore1_DEM)
dat_DEM_pth$lm_weights[dat_DEM_pth$nAttempts == 2 & dat_DEM_pth$cp == 1] <- 1/pscore2_DEM
dat_DEM_pth$lm_weights[dat_DEM_pth$nAttempts == 2 & dat_DEM_pth$cp == 0] <- 1/(1-pscore2_DEM)

did.out_boot_DEM_pth <- boot_lm(dat = dat_DEM_pth, formul = DV_diff ~ cp + factor(nAttempts), cluster_var = 'ID', strata = 'nAttempts')
l90_demo_pth <- quantile(did.out_boot_DEM_pth[,colnames(did.out_boot_DEM_pth)=='cp'], 0.05)
u90_demo_pth <- quantile(did.out_boot_DEM_pth[,colnames(did.out_boot_DEM_pth)=='cp'], 0.95)

did.out.DEM_pth <- lm(DV_diff ~ cp + factor(nAttempts),  dat_DEM_pth, weights = lm_weights)
coef_DEM_pth <- coef(did.out.DEM_pth)[names(coef(did.out.DEM_pth)) == 'cp']

coefs_pth <- c(coef_AR_pth, coef_DEM_pth)

l90_pth <- c(l90_autoc_pth, l90_demo_pth)
u90_pth <- c(u90_autoc_pth, u90_demo_pth)

pdf('FigureB1.pdf', width = 6, height = 4)
	par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(1.9,3.6,1.1,1.1))
	plot(x = c(1, 2), y = coefs_pth, xlim = c(0.5, 2.5), ylim = c(-0.1, 0.35), xaxt = "n", yaxt = "n", ylab = 'Effect of Coup Success', las = 2, xlab = '')
	axis(side = 1, at = c(1, 2), labels = c('Coups in Autocracies', 'Coups in Democracies'))
	axis(side = 2, at = c(0.4, 0.3, 0.2, 0.1, 0, -0.1, -0.2), labels = c(0.4, 0.3, 0.2, 0.1, 0, -0.1, -0.2), las = 3)
	abline(h = 0, lty = 2, lwd = 0.5)
	arrows(x0 = c(1, 2), y0 = l90_pth, y1 = u90_pth, angle = 90, code = 3)
dev.off()

#------------------------------
#Figure B2
#------------------------------
library(mice)
analyses <- list()
for (imp in seq(1, 10)){
	diff_dat_pth_t <- diff_dat_pth
	diff_dat_pth_t$latent.pt1 <- diff_dat_pth_t[[paste('posterior_latentmean', imp, sep = '')]]
	did.out.marg.ef.dat <- na.omit(diff_dat_pth_t[, c('DV_diff', 'cp', 'nAttempts', 'ccode', 'latent.pt1', 'Year', 'lm_weights')])
	did.out.marg.ef <- lm(DV_diff ~ cp*(latent.pt1) + factor(nAttempts), data = did.out.marg.ef.dat, weights = lm_weights)
	analyses[[imp]] <- did.out.marg.ef
}

pld <- pool(analyses)$pooled
rnms1 <- which(rownames(pld) == 'cp')
rnms2 <- which(rownames(pld) == 'cp:latent.pt1')
pld_cp_est <- pld[c(rnms1, rnms2), 'estimate']
pld_cp_var <- pld[c(rnms1, rnms2), 't']
pool_varVcov <- matrix(nrow = 2, ncol = 2, 0)
diag(pool_varVcov) <- (pld_cp_var)
draws <- mvrnorm(n = 100000, mu = pld_cp_est, Sigma = pool_varVcov)
colnames(draws) <- c('cp', 'cp:latent.pt1')

latentmean_levels <-  seq(-4, 4, 0.01)
frstTerm <- as.matrix(data.frame(rep(1, length(latentmean_levels)), latentmean_levels))
marginal_effects <- frstTerm %*% t(draws)
means <- apply(marginal_effects, 1, mean)
l90 <- apply(marginal_effects, 1, function(x) quantile(x, 0.05))
u90 <- apply(marginal_effects, 1, function(x) quantile(x, 0.95))
l95 <- apply(marginal_effects, 1, function(x) quantile(x, 0.025))
u95 <- apply(marginal_effects, 1, function(x) quantile(x, 0.975))

pdf('FigureB2.pdf', width = 16*0.5, height = 12*0.3)
	par(mgp=c(2, 0.5, 0), tcl=-0.4, mar=c(3.1,3.6,1.1,1.1))
	attach(mtcars)
	plot(x = latentmean_levels,
		y = means,
		type = 'l',
		xlim = c(-2, 3),
		ylim = c(-0.15, 0.33),
		xlab = 'Pre-Treatment Level of State Repression',
		ylab = 'Effect of Coup Success',
		xaxt = 'n',
		yaxt = 'n',
		lty = 1,
		lwd = 1
	)
	axis(1, at = -2:3, labels = -2:3)
	axis(2, at = seq(-0.2, 0.3, by = 0.1), labels = seq(-0.2, 0.3, by = 0.1), las = 3)
	polygon(c(latentmean_levels, rev(latentmean_levels)),
		c((l90), rev((u90))),
		border = FALSE,
		fillOddEven = TRUE,
		col = rgb(171/255, 171/255, 171/255, alpha= 0.55))
	polygon(c(latentmean_levels, rev(latentmean_levels)),
		c((l95), rev((u95))),
		border = FALSE,
		fillOddEven = TRUE,
		col = rgb(171/255, 171/255, 171/255, alpha= 0.25))
	abline(0, 0, lty =1, lwd = 0.5)
	rug(diff_dat$latentmean.pt1)
	abline(v = 0, lty = 2, col = 'grey')
	abline(v = 1, lty = 2, col = 'grey')
	abline(v = 2, lty = 2, col = 'grey')
	abline(v = 3, lty = 2, col = 'grey')
	abline(v = -1, lty = 2, col = 'grey')
	abline(v = -2, lty = 2, col = 'grey')
	abline(v = -3, lty = 2, col = 'grey')
dev.off()


#------------------------------
#Figure B4
#------------------------------
pth.keep <- !is.na(D_pth$cpAtL1)&!is.na(D_pth$latentmean)&!is.na(D_pth$latentmeanL2)
D_pth_tt <- D_pth[pth.keep,]
formul.B <- latentmean ~ cpAtL1*latent.pt1 + factor(Year) + trend*factor(ccode) + trendsqrd*factor(ccode) + trendcube*factor(ccode) # + trendfrth*factor(ccode) + shrYth1524
D_pth_tt$ccode <- as.factor(as.character(D_pth_tt$ccode))

coefs <- NULL
ses <- NULL
for (i in 1:10) {
	D_pth_tt$latent.pt1 <- D_pth_tt[[paste('latentmean_post_', i, '_L2', sep = '')]]
	D_pth_tt <- D_pth_tt[!is.na(D_pth_tt$latent.pt1), ]
	D_pth_tt$ccode <- factor(D_pth_tt$ccode)
	heter.out <- panel_reg(D_pth_tt, update(formul.B,  . ~ . + latent.pt1*cpAtL1), cus_weights = NULL )
	coefs <- rbind(coefs, heter.out$coefs[grepl('(^cp|^latent)', rownames(heter.out$coefs)),][,1])
	ses <- rbind(ses, heter.out$coefs[grepl('(^cp|^latent)', rownames(heter.out$coefs)),][,2])
	print(i)
}

c_pth <- as.vector(mi.meld(coefs, ses)$q.mi)
s_pth <- as.vector(mi.meld(coefs, ses)$se.mi)
vrNms <- colnames(mi.meld(coefs, ses)$q.mi)
names(c_pth) <- vrNms
names(s_pth) <- vrNms

rnms1 <- which(vrNms == 'cpAtL1')
rnms2 <- which(vrNms == 'cpAtL1:latent.pt1')
pld_cp_est <- c_pth[c(rnms1, rnms2)]
pld_cp_var <- s_pth[c(rnms1, rnms2)]
pool_varVcov <- matrix(nrow = 2, ncol = 2, 0)
diag(pool_varVcov) <- (pld_cp_var)^2
draws <- mvrnorm(n = 100000, mu = pld_cp_est, Sigma = pool_varVcov)
colnames(draws) <- c('cpAtL1', 'cpAtL1:latent.pt1')

latentmean_levels <-  seq(-4, 4, 0.01)
frstTerm <- as.matrix(data.frame(rep(1, length(latentmean_levels)), latentmean_levels))
marginal_effects <- frstTerm %*% t(draws[, grepl('^cp', colnames(draws))])
means <- apply(marginal_effects, 1, mean)
l90 <- apply(marginal_effects, 1, function(x) quantile(x, 0.05))
u90 <- apply(marginal_effects, 1, function(x) quantile(x, 0.95))
l95 <- apply(marginal_effects, 1, function(x) quantile(x, 0.025))
u95 <- apply(marginal_effects, 1, function(x) quantile(x, 0.975))

pdf('FigureB4.pdf', width = 16*0.5, height = 12*0.5)
	attach(mtcars)
	plot(x = latentmean_levels,
		y = means,
		type = 'l',
		#xlim = c(-2, 3),
		xlim = c(-2, 3),
		ylim = c(-0.15, 0.33),
		xlab = 'Pre-Treatment Level of State Repression',
		ylab = 'Marginal Effect of Coup Attempt on State Repression',
		xaxt = 'n',
		yaxt = 'n',
		lty = 1,
		lwd = 1
		)
	axis(1, at = -2:3, labels = -2:3)
	axis(2, at = seq(-0.2, 0.3, by = 0.1), labels = seq(-0.2, 0.3, by = 0.1), las = 1)
	polygon(c(latentmean_levels, rev(latentmean_levels)),
		c((l90), rev((u90))),
		border = FALSE,
		fillOddEven = TRUE,
		col = rgb(171/255, 171/255, 171/255, alpha= 0.55))
	polygon(c(latentmean_levels, rev(latentmean_levels)),
		c((l95), rev((u95))),
		border = FALSE,
		fillOddEven = TRUE,
		col = rgb(171/255, 171/255, 171/255, alpha= 0.25))
	abline(0, 0, lty =1, lwd = 0.5)
	abline(v = 0, lty = 2, col = 'grey')
	abline(v = 1, lty = 2, col = 'grey')
	abline(v = 2, lty = 2, col = 'grey')
	abline(v = 3, lty = 2, col = 'grey')
	abline(v = -1, lty = 2, col = 'grey')
	abline(v = -2, lty = 2, col = 'grey')
	abline(v = -3, lty = 2, col = 'grey')
dev.off()

